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We study decaying magnetohydrodynamics (MHD) turbulence stemming from the evolution of the Taylor-Green (TG) flow generalized 
"■^^ recently to MHD, with equal viscosity and magnetic resistivity and up to equivalent grid resolutions of 2048^ points. A pseudo-spectral 
I code is used in which the symmetries of the velocity and magnetic fields have been implemented, allowing for sizable savings in both 
^ computer time and usage of memory at a given Reynolds number. The flow is non-helical, and at initial time the Icinetic and magnetic 
^•j^ energies arc taken to be equal and concentrated in the large scales. After testing the validity of the method on grids of 512^ points, we 
analyze the data on the large grids up to Taylor Reynolds numbers of ~ 2200. We find that the global temporal evolution is accelerated 
in MHD, compared to the corresponding neutral fluid case. We also observe an interval of time when such conflgurations have quasi- 
O constant total dissipation, time during which statistical properties are determined after averaging over of the order of two turn-over times. 
* ^ A weak turbulence spectrum obtains which is also given in terms of its anisotropic components. Finally, we contrast the development 
j^.^ of small-scale eddies with two other initial conditions for the magnetic field and briefly discuss the structures that develop, and which 
^ display a complex array of current and vorticity sheets with clear rolling-up and folding. 

1 Introduction 

> 

Magnetic fields are present in many astrophysical and geophysical flows and are often important ener- 
2^ getically. Their origin - the dynamo problem, is ill-understood with many different behaviors observed in 
nature; for example, the solar magnetic field evolves in time in a somewhat regular fashion, leading to 
the possibility of a prediction of the amplitude and onset of the next cycle (Wang and Sheeley, 2006), 
whereas the terrestrial field has an erratic temporal behavior (Valet et al., 2005). The dynamo effect had 
^ escaped experimental verification until recently, due to the difficulty of reaching a sufficiently high magnetic 
Reynolds number for dynamo action to take place within a turbulent flow at low magnetic Prandtl number 
J> as occurs in liquid metals used in the laboratory; the study of such a held as produced in the experiment 
yields a wealth of information close to the threshold for dynamo action (Monchaux et al., 2007). However, 
;J] the high magnetic Reynolds number regime is still out of reach of the experimental approach using fluids; 
^ yet, due to its nonlincaritics, the full MHD problem at high kinetic and magnetic Reynolds number as 
encountered in the sun, the solar-terrestrial environment or the interstellar medium for example, is highly 
complex with many multi-scale interactions. 

It has been shown recently (Dar et al, 2001, Alexakis et al., 2005, Debliquy et al, 2005) that such 
interactions are noticeably more nonlocal than in the fluid case, involving widely separated scales. More- 
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2 MHD turbulence for Taylor-Green vortices 

over, several recent in situ observations in the magnetosphere and in the solar wind show the occurence 
of highly energetic events due to reconnection of magnetic field lines (Hasegawa et al., 2004, Sundkvist et 
al., 2005, Nykyry et al., 2006, Phan et al., 2006, Retino et al., 2007), as well as rotational discontinuities 
(Whang, 2004). It thus becomes important to be able to study in detail multi-scale interactions in MHD, 
a pre-requisite to which is to have sufficient scale separation in the fluid with high kinetic and magnetic 
Reynolds numbers. In three spatial dimensions, this represents a serious challenge from the numerical 
point of view, a challenge that will necessitate all the power that the world-wide petascale effort is going 
to offer, and more. There are many ways one can partially circumvent this difficulty, however; among them, 
modeling plays a central role. Enforcing numerically the symmetries that a given flow may have is one 
such way that wc shall employ in this paper. This allows for substantial savings in memory usage and in 
CPU time (although the CFL condition applies to the high Reynolds number that will be modeled this 
way). Using such highly symmetrical fields - that are extensions to MHD of the Taylor-Green (TG) flow 
studied first in the context of fluid turbulence (Brachet et al, 1983), we analyze the properties of MHD 
turbulence for several conflgurations (see Lee et al. (2008) for more details, and below) in a turbulent 
regime never explored before numerically in MHD in the incompressible case at such high resolution (see 
Vahala et al. (2008), and Kritsuk et al., 2009 for the case of very high resolution studies of compressible 
MHD turbulence) . The next section gives a description of the initial conditions; the temporal behavior of 
the insulating case (or IMTG hereafter) flow is given in §3; spectral behavior and structures are discussed 
in §4, and flnally Section §5 is the conclusion. 



2 Taylor-Green flows for MHD 

The MHD equations for an incompressible flow of constant unit density with a velocity v and magnetic 
induction h (in units of the Alfven velocity) read: 

dv 

— -\-v -Vv = -VP + j xh-\-u/S.v , (1) 



— = V X (v X 6) + r?A6 , (2) 

together with V-v = = V- 6;Pis the pressure and j = V x b the current density. In the absence of 
viscosity u and resistivity r), the total energy =< v^/2-|-6^/2 >= Ev + Em, the total cross-correlation 
He =< v b > and the total magnetic helicity Hm =< a • b > (where 6 = V X a is the magnetic potential) 
are conserved in three space dimensions. The MHD equations can take a more symmetric form using the 
Elsasser variables 2^ = v zt 6 with as invariants Hm and E± =< /2 >= Et ± Hq. Ideal flows in 
MHD have been studied numerically both in the two-dimensional case (see e.g. Prisch et al., 1983) and 
the three-dimensional case (Kerr and Brandenburg, 1999, Lee et al, 2008) including using adaptive mesh 
reflnement (Grauer and Marliani, 2000), but in this paper we concentrate on the dissipative case with 
a regular grid and using a fully parallelized pseudo-spectral code written specifically to implement the 
symmetries of the fiow (see below). Note that no uniform magnetic field is included in these computations, 
and there are no forcing terms either. 
The simplest Taylor-Green fiow can be written (Brachet et al., 1983): 

v{x, y, z) = vq [(sin x cos y cos z)ex — (cos x sin y cos z)ey\ ; (3) 

the velocity component in the third direction, equal to zero initially, will grow with time. As usual, the 
kinetic and magnetic Reynolds numbers are defined as 



Rv = vqL^Iv , Rm = VQL%jl'r] , 
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with iJyf^j, respectively the kinetic, magnetic and total integral scales L^j. = 2ttE^^ f[Ex{k)/k]dk with 
X = V, M, T and = j Ex{k)dk the kinetic, magnetic and total energy respectively with Et = Ey + Em- 
Similarly, one can evaluate Taylor Reynolds numbers RyMT based on the kinetic, magnetic and total 
Taylor scales defined respectively as: 

, ,j jEy,M,Tmk y/^ . 

these scales can be viewed as a measure of the curvature of the field lines. The flow is computed in a cubic 
box of length 27r with minimum and maximum wavenumbers kmin = 1 and kmax = -^/3 respectively, with 
N the number of points in each direction and using a standard 2/3 deliasing rule. The Taylor-Green vortex 
given above can be put in correspondence with the von Karnian flow between two counter-rotating cylinders 
as used in several laboratory experiments, including in the case of liquid metals (Bourgoin et al., 2002). The 
feasibility of dynamo action was shown numerically down to magnetic Prandtl number Pm = v/rj 10~^ 
using the TG flow in the kinematic regime in a combination of Direct Numerical Simulations (DNS) and 
modeling (Ponty et al, 2005), and a dynamo was recently obtained experimentally (Monchaux et al, 
2007, Bourgoin et al., 2007) using that configuration. The TG flow has also been used numerically to 
study sub-critical bifurcations in the dynamo regime (see e.g., Ponty et al., 2007), the hysteresis cycle 
being linked with changes in the velocity field associated with the Lorentz force. 

Generalizations of the Taylor-Green flow to MHD were presented in Lee et al. (2008) where the ideal 
{u = = r]) case was studied with initially the magnetic field taken as: 

b^. = feo cos(2;) sin(y) sin(z) (5) 
by = bo sin(x) cos(y) sin(2;) (6) 
bl = —2bo sin(a;) sin(y) cos(2;); (7) 

the velocity-magnetic field global correlation He = for this flow and remains so at all times because 
of the imposed symmetries, although it has been known for a long time that there may be strong local 
correlations corresponding to local alignment of the velocity and magnetic field (Passot et al., 1990). 
The magnetic induction is everywhere perpendicular to the walls (and therefore the current j = V x 6 
parallel to the walls) of the so-called impermeable box defined as [0, vr]^; the impermeable box thus appears 
to be insulating, and this fiow is henceforth being named Insulating Magnetic Taylor-Green fiow (IMTG 
hereafter). 

An alternate initial conditions for the magnetic field, hereafter, is insulating as well: 

b^ = 6o cos(2x) sin(2j/) sin(2z) , by = -6o sin(2x) cos(2y) sin(2z) , bf = 0, 

with again zero total cross helicity. Finally, one can also construct a set of initial conditions, labeled 
"conducting" (the current being perpendicular to the impermeable box): 

6^ = b^ sin(2x) cos(2y) cos(22:) , b^ = b^ cos(2a;) sin(2y) cos(22;) , 6^ = -26^ cos(2a;) cos(2y) sin(22;). 

In this latter configuration {Be hereafter), H'-' is non-zero but very weak (less than 4% at its maximum, 
a dimensionless measure of correlation, relative to the total energy), whereas the magnetic helicity is zero 
in all three configurations. The parameters vq and 6o are chosen so that at initial time, the kinetic and 
magnetic energies are equal with Et = 0.25. 

Within the periodic cube of length 2it, mirror symmetries about the planes x = 0, x = ir, y = 0, y = 
IT, z = 0, and z = n are present and will be enforced numerically, together with rotational symmetries of 
angle Nn about the axes (x, y, z) = (f , y, f ) and (x, |, |) and of angle Ntt/2 about the axis (|, |, z), for 
N e Z (see Brachet et al., 1983, for more details). Note that the flow has several nulls (three components 
equal to zero): the central point to the impermeable box Ni = {it/2, 7r/2, 7r/2), and the planes Pi = {x, 0, 0), 
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Figure 1. Left: Flux-law for the z+ = v + b variable (see eq. 10), compensated by i: the exact law appears as a horizontal line. 
Right: Probability distribution function of cos[v • b]. Data for both figures are taken near the peak of dissipation for 
computations on grids of 512* points, with the full run drawn with a solid line and the symmetric run with diamonds. 

P2 = (0, y, 0) and P3 = (0, 0, z). Partial nulls (one or two components equal to zero) also occur; they may 
be as well the sites of formation of strong current sheets and thus of reconnection events. 



3 Temporal behavior of the insulating TG flow (IMTG) 
3.1 Is a symmetric run faithful to the full dynamics? 

The numerical method we employ in this paper is based on the assumption that the implementation of 
the TG symmetries will not alter the overall evolution of the flow. We show in Figs. 1 and 2 that this is 
indeed the case, at least at the grid resolution employed for this test, using 512^ points. The data is taken 
at peak of dissipation, once the turbulence has fully developed. 

The first test we performed is to verify the exact law that can be written in MHD in terms of third-order 
structure functions (Politano and Pouquet, 1998); in dimension three they read: 

( SvL {5vif ) + {5vL {5hif ) - 2( 5hL Svi dbi) = -^e^ r , (8) 

-( 5bL (Sbif )-{5hL {6vif ) + 2( 5vL Svi Sh ) = r ; (9) 

these laws can be written in a more compact form, using the Elsasser variables z='= = vitb, with associated 
energies 

{5zl{v) [5zf{v)]^) = -^6± r , (10) 

where = c'^ it e'~' arc the energy transfer rates of E^^ and and e*^ the rates for total energy and 
velocity-magnetic field correlation. These laws, under the assumptions of incompressibility, homogeneity, 
isotropy, stationarity and large Reynolds number are nothing more than the expression of the invariance of 
quadratic functionals of the fields under the MHD dynamical equations; as usual, Sfi{r) = /j(x + r) — /j(x) 
is the difference for the i*'*-component /j of the field /, and is its longitudinal component, i.e. the vector 
field projected on the direction r along which the difference is taken. The laws in terms of the Elsasser 
fields show that the two non-helical invariants, E^ or equivalently Et and H'-'\ are coupled; this leads to 
a double direct cascade towards small scale. In terms of the velocity and magnetic field, it shows that the 
conservation of H'-^ is on equal par with that of total energy Et and that a priori two time scales can be 
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Figure 2. Visualization, using the VAPOR software (Clyne et ah, 2007) of the current density, on grids of 512^ points. Full DNS 

(right) and run implementing symmetries (left). 

expected to play a role in MHD dynamics, associated with the two invariants which are known to provide 
a partition of phase space (together with magnetic helicity), as shown in Stribling and Matthaeus (1991). 

We see in Fig. 1 (left) that the flux law is verified on a small interval of scales, representing the inertial 
range at that grid resolution of 512^ points, and that the full and symmetric runs give identical results. 
Similarly, when comparing the Probability Distribution Functions (PDFs) of the alignment between the 
velocity and the magnetic field, as measured by the cosine of the angle between the two vectors, again no 
distinction can be made between the full and symmetric run: they both have a central peak (corresponding 
to orthogonality of v and b) due to a persistence of initial conditions (not shown), and a dynamical tendency 
toward alignment, as observed in many MHD flows (see for a recent discussion, Matthaeus et al., 2008, 
Servidio et al., 2008). Finally, a three-dimensional snapshot of the current density (Fig. 2) for the whole 
computational box at peak of dissipation shows that the same type of structures appear in configuration 
space, with a predominance of sheets, some with strong curvature. 

3.2 Energetics of the IMTG flow 

We now describe the overall energetics of the IMTG initial condition given in Eqs. (5,6,7) in the presence of 
dissipation at high Reynolds number. At first, the ideal behavior analyzed in Lee et al. (2008) is recovered 
with an exponential decay of small scales corresponding to the formation of current and vorticity sheets, 
and a quasi conservation of total energy, even though an exchange between the kinetic and magnetic energy 
Ey and Em is occurring as can be seen in Fig. 3; this exchange can be related to an Alfvenic effect based 
on the large-scale magnetic field (recall that there is no imposed uniform field in these runs). The energy 
ratio Em / Ey is larger than unity at all times for this flow; a similar evolution obtains for the alternate flow 
Ba whereas, in the case of the conducting flow Be, the ratio is smaller than unity after a short transient 
(see Table 1 for more details). When displaying the temporal evolution in log-log coordinates for the total 
energy of the IMTG flow, an approximate power-law decay can be observed for a short time, with an index 
clearly smaller than what a Kolmogorov analysis would give: the energy decay in the Navier-Stokes case 
is expected to vary as {t — t*)"^'^/'' with t=K a characteristic time of decay of the energy, of order unity here 
given our initial conditions; however, in the case of the TG fluid flow, it has been found (Brachet et al., 
1983) that energy decay follows a {t — t^)'"^ power law, a law that can be recovered phenomenologically 
using an argument based on the fact that the integral scale cannot grow in the TG flow used here because 
the initial conditions are centered on the largest available scale (Cichowlas et al., 2005). 

Slower temporal behavior in MHD, compared to the pure fluid case, has been observed by a number 
of authors (Hossain et al, 1995, Kinney et al, 1995, Galtier et al, 1997, 1999, Bigot et al, 2008). It 
can be attributed to the slowing-down of energy transfer to small scales due to the interactions between 
waves and turbulent eddies as modeled originally by Iroshnikov (1963) and Kraichnan (1965) under the 
simplifying assumption of global isotropy (see Goldreich and Sridhar, 1997, and Ng and Bhattacharjee, 



June 7, 2009 20:5 Geophysical and Astrophysical Fluid Dynamics gafd'MHD'v5 

6 MHD turbulence for Taylor- Green vortices 

0.3 



0.2 - 

LjJ 

0.1 - 



0.0 

2 4 6 1 2345678 

t t 

Figure 3. Left: Total (solid), kinetic (dot) and magnetic (dash) energies as a function of time for the IMTG flow; equivalent grid of 
2048^ points, Taylor Reynolds number of Rx ~ 2200; note the excess of magnetic energy at all times, and the slow decay of the energy 
compared to the fluid case. Right: In log-log coordinates, the total energy decay appears to follow a weaJj power law. 

1997, for a straightforward extension of that phenomenology to the anisotropic case). This slowing-down 
can be understood when writing the MHD equations in terms of the Elsasser variables z='= = v ± b: the 
nonlinearities involve the products ■ Vz^ but there are no self interactions {z'^z'^ or z~z~). In MHD, 
following the same phenomenology as Kolmogorov but taking into account the effect of Alfven waves, 
the temporal decay of energy can be shown to become, in three dimensions ~ (t — t^,)^^^^ under such 
assumptions. Note that a further slowing-down of energy decay may stem from the presence of a strong 
large-scale field, as shown recently in the context of anisotropic MHD (Bigot et al, 2008), with now a 
decay ~ — t*)~^/^. However, the slowing-down observed here, with E{t) ~ (t* — t)~^'^, is even more 
important. The origin of this behavior is not known. It may be related to the strong (relative) increase 
of magnetic energy, with Em/ Ey reaching a maximum ~ 4 at t = 2.5 and decreasing steadily thereafter 
to a value close to 1.5 at the end of the run; this may entail a faster decay at later times that could 
only be observed at higher Reynolds numbers, although this peak in the magnetic to kinetic energy ratio 
augments monotonically with Reynolds number, from Em/Ev ~ 2.5 for R?" ~ 120 (see also Table 1). 
This slow decay could also be due to the local alignment between the velocity and magnetic field although 
the total correlation remains weak in relative terms (less than 4%). Tendency toward alignment of v and 
b is obtained in many observational and numerical MHD flows as mentioned earlier (Matthaeus et al., 
2008, and references therein); it weakens the nonlinear terms responsible for the decay of energy at high 
Reynolds number; indeed, slow decay can be observed in the presence of strong and global correlations 
(Galtier et al., 1999). Another feature specific to the IMTG flow is that, at t = 0, the magnetic field and 
the vorticity are identical (and anti-parallel). Such is not the case for the other two flows studied in this 
paper, for which a faster decay obtains with the total energy varying as {t — t*)^^"^, although still slower 
than for the pure fluid TO case. Obviously higher resolutions would be needed to clarify this point for the 
IMTG flow since it would allow for a larger temporal range in the self-similar decay of the energy before 
the energy and thus the Reynolds numbers drop too significantly. 

When examining the temporal evolution of the maxima of the current (dash line) and the vorticity (dot- 
ted line), given in Fig. 4 (left), we observe that the first initial exponential phase is ideal and corresponds to 
thinning of current and vorticity sheets due to large-scale shear; this phase is interrupted rather abruptly 
at t 2.5 by another phenomenon that now evolves more rapidly and corresponds to the quasi-collision 
of two current sheets, pushed together by the contribution of the Lorentz force to the magnetic pressure 
on either sides of the sheets (see also Fig. 9 in the next Section). This is still in the ideal non-dissipative 
phase and gives rise to a quasi rotational discontinuity (Lee et al., 2008) which is then arrested, around 
t 2.65, by dissipation setting in at the smallest scales. The dissipative phase in the sup norms can be 
seen as somewhat chaotic but a characteristic feature emerges, that is that the ratio ^maxl^max remains 
approximately constant, including at the latest time of the computation when roughly 28% of the energy 
has already been lost to dissipation. This is somewhat reminiscent of the constancy of the cancellation 
exponent k in a turbulent flow (Graham et al., 2005), where k measures by how much the change in sign 
of (say) the magnetic field at a given scale varies with scale; it can be attributed to the fact that in the 
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Figure 4. Left: Maxima of the current (dash line) and vorticity (dotted line) as a function of time, with a slight excess of the former. 

Right: Temporal evolution of dissipation (kinetic, dotted line, magnetic, dash line, and total, solid line) as a function of time. 
Note the plateau first reax;hed around t « 3.5 and lasting for a couple of eddy turn-over times, and subsequent peaks likely 

associated with reconnection events of current sheets. 



self-similar decay of energy, as long as the Reynolds number remains sufficiently large, the complexity of 
the flow (as measured for example by k) remains approximately the same even though the energy itself 
does decay, at an algebraic slow rate. 

In Fig. 4 (right) is shown as a function of time the total generalized enstrophy Jl^ =< > + < > 
(solid line), proportional to energy dissipation since u = ij, as well as its kinetic (dots) and magnetic (dash) 
components. After an initial exponential growth, the IMTG flow displays a sizable plateau during which the 
dissipation of energy remains quasi constant and thus during which temporal averages can be performed 
in order to study with greater accuracy the statistics of the flow such as its spectral behavior (see below). 
Again, the small scales are dominated by the magnetic structures; also note the presence of later maxima 
at t 5 and t k;7, probably a sign of renewed reconnection events of magnetic field lines; the second peak 
around t = 5 appears at this resolution but is barely perceptible at lower Reynolds numbers. 




0246 02468 

t t 

Figure 5. Left: Temporal variation of the magnetic to kinetic ratios of the first three moments of the basic fields (see Eq. 11) with 
n = 0, 1, 2 represented with a solid, dash and dotted line respectively. All display an excess of magnetic excitation, with an earlier peak 
for those moments involving derivatives and corresponding to an abrupt acceleration in the development of small scales in that flow. 
Right: Temporal variation as a function of Reynolds number (and thus grid resolution) of total dissipation for the IMTG flow, from a 
grid of 128"^ points (long dash) to 2048"^ (dots) in lin-log coordinates; the inset in log-log coordinates gives the first peak of dissipation 
as a function of Ry; note the tendency for it to slowly decrease with Reynolds number for the IMTG flow. 



This is corroborated by the examination of the ratio of magnetic to kinetic excitation for several moments 
of the fields displayed in Fig. 5 (left); these ratios are defined as: 

" Jk^^Ev{k)dk 

and can be related to the respective Sobolev norms. There is a systematic overshoot of magnetic excitation, 
a well-known feature of turbulent flows for example often observed in the solar wind for the energy ratio 



(11) 
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Ro; it may possibly be due to the fact that turbulent magnetic resistivity is known to be less efficient 
than turbulent viscosity (Pouquet et al., 1976), a phenomenon modeled through the use of second-order 
closures of turbulence. Small-scale structures are dominated as well by magnetic excitation, to a somewhat 
lesser extent as n increases and with an earlier peak at t ~ 2.5 when the first strong current structure has 
formed (see Lee et al., 2008, and Fig. 9). Furthermore, Rn 1 as n increases (and also as time elapses); 
this could be interpreted as being due to the fact that, as we increase the weight of the small scales, these 
global averages get to be dominated by a few localized events in space that tend to be similar in their 
physical structure. 

The variation of total dissipation with time of the IMTG flow is given in Fig. 5 (right) in lin-log 
coordinates for several Reynolds numbers. We observe that the initial phase shows less dissipation when 
the viscosity is reduced (from 2 x 10^^ (long dash) to 6.25 x 10~^ (solid line), with runs on grids of 128'^ 
points to 2048^ points, augmenting by a factor of 2® and with two different runs at the largest resolution. 
The inset gives, in log-log coordinates, the variation with Reynolds number of the first temporal maximum 
of dissipation at the end of the ideal phase; for this flow, there seems to be a steady (power-law) decrease 
of dissipation, although constant energy dissipation is observed in MHD turbulence with the Ba,c flows, 
as well as in a full numerical simulations using a Beltrami flow (Mininni and Pouquet, 2007). The question 
of finite dissipation in the limit of infinite Reynolds number in MHD will thus need further study at 
higher Reynolds number; indeed, one difference in behavior between the IMTG and Ba,c flows is the afore 
mentioned observation of a very slow decay of energy for IMTG when contrasted to the other flows. 



4 Spectra cind structures for the IMTG flow 
4.1 Energy spectra 



0.01 





Figure 6. Total energy spectra compensated by fc™, with m = 2, 5/3 or 3/2 (see labels) for the IMTG flow. Left: instantaneous total 
energy near the maximum of dissipation. Right: total energy spectra averaged in the plateau of dissipation, t £ [3.5, 5]. Note that in the 
figure at right, the spectra are also averaged between adjacent shells in order to get rid of the even-odd oscillations due to the structure 
of the Taylor-Green configuration. A clear spectrum emerges from this analysis for the total energy. 



We show in Fig. 6 (left) the total isotropic energy spectrum, compensated by (top curve), 
(middle curve) and fe^/^ (bottom curve), computed at i = 4 close to the first peak of dissipation; these 
compensations correspond respectively to a weak turbulence spectrum (WT hereafter, Galtier et al., 2000, 
2002, 2005), to a Kolmogorov (1941) spectrum (hereafter, K41) and to an Iroshnikov (1963) - Kraichnan 
(1965) spectrum (hereafter, IK). Issues relating to anisotropy are discussed further below. The spectrum 
is instantaneous and it is computed for the run on a grid corresponding to 2048^ points. Although even at 
such a resolution, it may be hard to discern spectral laws, the WT spectrum seems to be fulfilled better. Fig. 
6 (right) is the same as Fig. 6 (left) but now the data is averaged over the plateau of energy dissipation 
discussed previously and over adjacent shells; a k~'^ law emerges clearly when temporal averaging is 
employed. Note that the kinetic and magnetic integral scales for the IMTG fiow 16 (see Table 1) are 
respectively 1.61 and 1.95, the Taylor scales are 0.32 and 0.42 and the Kolmogorov dissipation scales 
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0.03 and 0.02. These values obtain from a temporal averaging between t=4.5 and t=5, over 21 snapshots; 
we thus conclude that the flow is well resolved since the dissipation scale is measurably larger than the 
grid spacing, and that the inertial index for the total energy is a nonlinear effect, followed by a sizable 
dissipation range. Further, we note that the Alfven time is 3.3 times shorter than the eddy turn-over time 
at the integral scale, and is almost 14 times smaller at the Taylor scale, consistent with the fact that the 
spectra show no K41 behavior, the transfer time of energy through nonlinear mode coupling being affected 
by Alfven wave propagation. 



1 .oo 



O-O 1 

o.o 1 o. 1 o 1 .oo 

Figure 7. Perpendicular (solid line) and parallel (dotted line) components of the second-order structure function of the magnetic field, 
compensated both by I; ± and || refer here, in the absence of an imposed magnetic field, to the direction of a locally-averaged magnetic 
field, in a sphere of diameter the magnetic integral scale. Note that the perpendicular structure function follows a law 52, x ~ I over 
roughly a decade in scales, corresponding to a weaJs turbulence regime (Galtier et al., 2000). Similar results obtain for the velocity. 

In the case of weak turbulence, it has been known for a long time that the flow becomes anisotropic 
under the influence of a large-scale (uniform) magnetic fleld Bq; a phenomenological description of the 
influence a large scale magnetic fleld can have on the anisotropic dependence of spectra (Goldreich and 
Sridhar, 1997, Ng and Bhattacharjee, 1997) yields, by a simple generalization of the IK reasoning, to a 
fcj^ spectrum, as conflrmed by theoretical developments (Galtier et al., 2000). 

ADD HERE: Weak turbulence has been observed in the magnetosphere of Jupiter (Saur et al., 2002) 
and more recently in a large numerical simulation (Mininni and Pouquet, 2007). Numerous studies of 
two-dimensional MHD turbulence, which can be regarded as a simple model of MHD turbulence in the 
presence of a strong magnetic fleld, have been performed, including in its intermittency properties (see e.g. 
Sorriso et al., 2000); the extension of this work to the study of non gaussianity of weak MHD turbulence 
is left for future work (Lee et al, 2009). 

In the absence of a uniform fleld, the local mean fleld at large scale can play an equivalent role as far as 
the small scales are concerned, provided there is sufficient scale separation between the two, leading to a 
strong energetic difference between the largest and smallest scales, by several order of magnitudes. High- 
resolution runs provide, to some extent, such a scale separation and we now analyze the flow in these terms, 
following the procedure developed previously (Mininni and Pouquet, 2007) and applying it to the magnetic 
energy spectrum: speciflcally, one averages the magnetic fleld in a sphere of diameter the magnetic integral 
scale, and deflnes, locally, the perpendicular and parallel dependence of structure functions with respect 
to that locally averaged fleld. The result can be found in Fig. 7 with the second-order structure function of 
the magnetic fleld displayed in terms of its perpendicular (solid line) and parallel (dotted line) components, 
and both compensated by i (i.e., corresponding to a k7\ spectrum); these functions are plotted at t = 5, 

-'-'II 

at the second peak of dissipation (see Fig. 4). A convincing weak turbulence scaling again appears in the 
large scales for the energy, in terms of l±, starting roughly at the magnetic Taylor scale (Am ~ 0.42), 
whereas the smallest scales follow a regular behavior as given by a Taylor expansion, again a testimony of 
the fact that the flow is well resolved. The structure function in terms of its variation with , regular again 
in the smallest scales, has a more chaotic behavior in the large scales and does not seem to follow a clear 
power law. It should be noted that the theory for weak MHD turbulence does not give any dependence 
on which thus may arise from higher-order developments. Considering the phenomenology proposed 
by Goldreich and Sridhar (1995), of a scale-independent equilibrium between the timescales for nonlinear 
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eddy turn-over time and the Alfven time, is a further step in the description of such flows, a point that 
win be studied in the future (Lee et al., 2009). 

Several other features in these plots are noticeable. First of all, the magnetic field is dominated ener- 
getically at all scales by its parallel component, a feature likely dependent on both the initial conditions 
and the intrinsic dynamics of the flow. Furthermore, the ratio of the parallel to perpendicular structure 
function is somewhat constant at large scale, and thus it appears that there is similar transfer at large 
scales (scales larger than the magnetic Taylor scale) in the perpendicular and parallel components of the 
energy spectrum; this confirms the finding of Mininni and Pouquet (2007) where, due to initial conditions 
that are strongly isotropic, the large scale spectrum remains isotropic down to the magnetic Taylor scale, 
and only becomes compatible with weak turbulence at scales smaller than Am- 

4.2 Small-scale structures 

The abrupt transition in the scaling behavior of the second-order structure function of the magnetic field 
just discussed takes place close to the dissipation scale; this can be viewed as being indicative of the 
formation of sharp structures in configuration space. We give in Fig. 8 a perspective volume rendering, 
zooming on different structures of the IMTG flow at time t = 5. The opacity for the plots is such that 
only features above a few r.m.s. value are displayed. On the left is a current sheet which has formed a roll; 
the blue and purple lines indicate magnetic field lines which display the formation of a cusp-like feature 
close to the core of the rolling current sheet; this roll is roughly 7r/20 across, which seems typical when 
one examines other current sheets in the computational box, with some variation in transverse size and 
with a length of the order of the magnetic integral scale. In Fig. 8 (middle), the vorticity is plotted at the 
same time and at the same location in space. This confirms earlier findings (Mininni and Pouquet, 2007) 
of strong correlation between the velocity and magnetic field in the small scales (as characterized by the 
behavior of vorticity and current), even though the global correlation remains small (of the order of a few 
% as mentioned before); this implies that, with weak nonlinear terms in such structures, they will survive 
for a time of the order of an eddy turnover time, similar to the case of a fluid within which the vortex 
filaments display a strong correlation between velocity and vorticity. Note that the curl of the Elsasser 
variables, lo^ = lo ±j, are also co-located, with a dominance of \Lii~ \ by 36%. 



Figure 8. Perspective volume rendering of a zoom on the intensity of the current (left), the vorticity (center) and another sub- volume 
for the current (right), all taken near peak of dissipation for the IMTG flow on a grid of 2048^ points, with a Taylor Reynolds number 
of ~ 2200. Note the curved vorticity and current sheets that are highly correlated spatially, and the sharp roUing-up of field lines in the 
current sheet. At right, it is shown that the Kelvin-Helmoltz instability of the current sheet (in a different sub-volume of the flow) can 

occur with different wavelengths. 

On the right of Fig. 8, a zoom in the current is shown, at a different location, with a small-scale wave-like 
structure in the sheet, this time with a characteristic scale rather like 7r/50, still well resolved numerically. 
These results confirm the instability, at high Reynolds number, of vorticity and current sheets that form 
into rolls; it also shows, possibly for the first time in direct numerical simulations of MHD, a small-scale 
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Figure 9. Temporal evolution of the logarithmic decrement 5 (see Eq. 12) in units of maximal wavenumber kmax for runs 16 (solid 
line), A6 (dash line) and C6 (dotted line), showing that each run is well resolved. See Table 1 for the nomenclature of the runs. The plot 
shows an abrupt event for t Ri 2.5 for the 16 flow (earlier for the magnetic skewness), linked to the quasi-coUision of two current sheets. 



undulating structure likely linked to a Kelvin-Helmoltz instability, but at a slightly different scale than 
the rolls on the left. Such Kelvin-Helmoltz instabilities have been observed in the magnetosphere (see e.g. 
Hasegawa et al, 2004, Nykyri et al., 2006) in a much more complex physical settings than what is being 
studied here, with compressibility, boundary and plasma kinetic effects in particular. The occurence of 
a Kelvin-Helmoltz instability of current sheets has been widely studied in the presence of an imposed 
shear (see e.g. Knoll, 2002, 2003), but it has being obtained only recently numerically in three dimensions 
in a fully turbulent incompressible MHD flow. Recent results using the Piecewise Parabolic Method for 
supersonic MHD turbulence also show roUcd-up current sheets, in a series of studies performed in the 
context of modeling the interstellar medium (Kritsuk and Padoan, 2009). 

One of the issue concerning the development of structures in a turbulent flow is to what degree the 
relevant vector fields arc aligned; indeed, the nonlinearities in the MHD equations can be written in term 
of the Lamb vector u x u;, the Lorentz force j x b and Ohm's law v x b. It can be shown that the primitive 
equations lead to an enhancement of the alignment between such vectors when they involve invariants such 
as the kinetic helicity (in the pure fluid case) or the cross-correlation between the velocity and magnetic 
field; such a local (as opposed to global) spatial alignment between the velocity and magnetic field is clearly 
observed both in the Solar Wind (using Ulysses data) and in DNS (Matthaeus et al., 2008, Servidio et 
al., 2008). In the case of the flows studied in the present paper, alignment obtains again, but with also 
a peak in the PDF at orthogonality of vectors, a feature directly linked to the chosen initial conditions. 
However, a more detailed analysis reveals that the Lorentz force retains its full strength in all our TG runs 
(particularly so for the Ba flows), in sharp contrast with the results of Servidio et al, 2008; moreover, the 
IMTG and Ba flows both show a tendency toward Alfvenic alignment (v || b) while retaining a peak at 
orthogonality due to initial conditions (sec Fig. 2, left), whereas the Be flow displays a rather flat PDF, 
in contrast to the other two flows. Such alignment properties are not likely due to the enforcement of 
symmetries: Fig. 2 (left) gives the same data for the symmetric run (solid line) and the full DNS (circles) 
on a resolution of 512^ grid points, with no discernable discrepancy. 



4.3 Comparative behavior of the three initial conditions 

We have emphasized the properties of the IMTG flow until now. In this section, we want to examine briefly 
the development of small-scale turbulence for the other two types of initial conditions that we have written 
in §2 for a symmetric Taylor-Green MHD code; Table 1 summarizes several of the features of such flows, 
for the lowest to the highest Reynolds numbers. 

One possible way to monitor the development of small scales is to examine the temporal development 
of the logarithmic decrement. The isotropic total energy spectrum ET{k,t) defined in the usual manner 
by averaging in spherical shells of width Ak = 1 can be fitted with the following functional form (Sulem 
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Table 1. Parameters of the runs described in this paper. The type of run (IMTG, Ba or Be is indicated respectively as I, A or C, with run "1" 
or "2" for computations at a grid resolution of 128^ points, and runs labeled "6" on equivalent grids of 2048^ points. All flows have equal kinetic 
and magnetic energy Ey — Em — 0.125 at t = 0; viscosity for run II is 2 x 10~^, for runs A2 and C2 , it is 10~^, and for runs 16, A6 and C6, 
it is equal to 6.25 X 10~^, with in all cases v — rj. Rv and R'^j are respectively the Reynolds number at initial time and the magneticTaylor 
Reynolds number at the final time of the computation; 6 is the logarithmic decrement whose expression is given in Eq. (12) and evaluated on 
the total energy spectrum; kmax is the maximum wavenumber, with N the linear number of points in the grid (here, 128 or 2048). The ratio 
of magnetic to kinetic energy Em/Ev is the maximum value reached by this ratio as it evolves with time (excluding t — 0). The magnetic to 
kinetic modal ratio = EM{kx)/ (l^x) is evaluated as an average over of the order of two turn-over times after the peak of dissipation, with 
X = 2,1 and A for wavenumbers fe = 2, kj^ and respectively, k^^^ being the magnetic integral and Taylor wavenumber (see Eq. 4 for the 
associated length scales). Similarly, Qm /^v =< > / < o?^ > is the maximum value attained by the ratio of magnetic to kinetic dissipation 
over time (recall that f = tj). Finally, S'^'^^ and 5^*"''*' are the skewness evaluated on the velocity field at the time of maximum dissipation and 
the final time respectively. 

* Note that in the case of the Be runs, the ratio of total magnetic to kinetic energies are given at their minimum value in time (this global ratio 
never goes above unity as time evolves). ** For the Be enstrophies, the values are those at the first maximum after the initial time. 



Run 


Rv 


7?A 




Em/Ev 




Pf 


pf 




gmax 


nfinal 


11 
16 


800 

2.3 X lO'' 


110 
2200 


1.9 
2.2 


2.7 
6 


19 
3 


7 

28 


2 

18 


5.7 
6.1 


1.5 
2.8 


1.15 

1.15 


A2 
A6 


1200 
2 X 10^ 


160 
900 


2.1 

3 


1.5 
1.9 


0.3 
0.4 


1.9 
2.2 


1.8 
1.6 


2.2 
4.4 


2.2 
2.2 


0.6 
0.8 


C2 
C6 


1200 
2 X 10^ 


100 
900 


0.8 
2.8 


0.35 * 
0.65 * 


0.09 
0.06 


1.2 
2.2 


1.7 
1.5 


2 * * 
2 


2.1 
2.8 


0.55 
1.6 



et al, 1983): 

ET{k, t) = - {[v{k, t)f + [b{k, t)]'^) d^k = C(f)fc-"We-2^W'= . (12) 

2 7|fc| 

The so-called logarithmic decrement S is a measure of the smaUest scale reached by the flow at a given time 
and can be compared to the grid resolution 3/N (taking ahasing into account); this criterion thus allows for 
measuring the accuracy of the computation at any given time, comparing 5 and 3/N. The analysis of the 
temporal evolution of 6 was performed in the ideal (non dissipative) case both in two dimensions (Frisch et 
al., 1983) and in three dimensions (Lee et al., 2008). In the latter case, it is found to decay exponentially 
{S{t) = 6oe~^/'^*), with two different values of Tj corresponding to two different ideal instabilities: first, a 
classical thinning of current and vorticity sheets occurs, due to the large-scale effect of shear; this phase 
is followed by a faster near-collision of two near-by current sheets, due to a favorable magnetic pressure 
gradient associated with the Lorentz force. The computation in the ideal case has to be stopped once 
S ~ but in the dissipative case, the run can be continued for long times provided the viscous and 

resistive terms are strong enough to arrest the process of small-scale formation. 

In this context, we now examine the temporal behavior of the logarithmic decrement in the dissipative 
case (see Fig. 9). We observe that at the same time as for the ideal case, an abrupt change occurs for 
the IMTG flow (solid line) , followed by a saturation of the decreasing of 5 when dissipation sets in, at a 
wavenumber wich, in the fluid case, varies as Ry'^- Only the IMTG flow displays a substantial acceleration 
of the development of small scales, also discernible in the evolution of the normalized third order moment 
of the velocity and magnetic field derivatives, but the three flows are well resolved at all times since 
Skmax > 2. The development of small scales appears more rapid for the i?A,c flows, since they reach their 
plateau in S at an earlier time (~ 2) than for the IMTG flow. Similar behavior can be observed on the 
temporal evolution of the kinetic and magnetic skewness (normalized third-order moment of the velocity 
and magnetic field), as seen in Fig. 10. Note that these moments systematically reach higher values than 
in the pure fluid case, and that one finds for all three flows higher values of the velocity skewness compared 
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Figure 10. Left: Temporal evolution of the skewness, i.e. the normalized third-order moment of the velocity field, all flows computed 
on an equivalent grid of 2048^ points, for runs 16 (solid line), A6 (dash line) and C6 (dotted line) (see Table 1). Right: Skewness built 
on the magnetic field for the same three runs. Both plots show an abrupt event for t Ri 2.5 for the IMTG flow (earlier for the magnetic 
skewness) , linked to the near-collision of two current sheets leading to a quasi rotational discontinuity. 



to its magnetic counterpart. 

The three types of runs studied here have almost identical initial conditions, from an ideal point of view: 
same energy, same equipartition at t = between kinetic and magnetic energy, same zero magnetic helicity, 
same weak correlation between the velocity and the magnetic field (0 to 4% in relative terms). They also 
have the same type of structures in the current and vorticity, with sheets and rolls. However, they differ 
in their dynamics. For example, for the IMTG runs, the magnetic to kinetic energy and enstrophy ratios 
are all larger than unity, particularly so at A; = 2 (and the length scales built on the velocity are smaller 
than those build on the magnetic field). For the Ba runs, the kinetic and magnetic variables are close to 
being in balance whereas in the Be case, it is now the kinetic quantities that dominate. The strongest 
differences between the three flows is the dynamics of the low-k modes, whereas, for all runs, there is an 
excess of magnetic energy for k > kf^ , with kf^ the magnetic integral wavenumber (see Table 1). These 
issues will be examined further in the future (Lee et al, 2009). 



5 Discussion and conclusion 

In this paper, we have implemented numerically the symmetries of the Taylor-Green flow generalized to 
MHD and have thus been able to study the dynamics of MHD turbulence at higher Reynolds number 
than what can be reached in a full DNS computation. We checked that no spurious error develops when 
using such codes, and we have explored the properties of three types of initial conditions, stressing the 
evolution in one particular case, the IMTG flow. We show that this flow develops an energy spectrum, 
Exik) ~ which is in agreement with the prediction of weak MHD turbulence, the spectrum being 
obtained as well for Exik^^ k'J^ . The detailed analysis of the other two sets of initial conditions that 
follow the symmetries of the Taylor-Green flow are left for future work (Lee et al., 2009). 

In the study of MHD turbulence and of the dynamo problem, the Taylor-Green flow has proved quite 
useful. Other methods to progress in our understanding of MHD turbulence is to resort to modeling. 
Andrew Soward (1972) pioneered this in developing an Eulerian-Lagrangian approach to the kinematic 
dynamo (see also Soward and Roberts, 2008, for a recent analysis). Other models can be derived and 
used (see e.g., Kraichnan and Nagarajan, 1967, Pouquet et al., 1976, Yoshizawa 1990, Holm 2002a, 2002b, 
Graham et al., 2009), including numerical (see e.g. Meneguzzi et al., 1996). Using both the implementation 
of symmetries and such models may prove a fruitful approach to explore parameter space (Pouquet et al., 
2009). For example, it has been known for a long time that the amount of correlation between the velocity 
and the magnetic field modifies the spectral indices of the Elsasser variables (and hence of the other energy 
spectra), with a steeper slope for the dominant mode (see Grappin et al., 1983, in the context of isotropic 
turbulent closures, Galtier et al., 2000, for a weak turbulence anisotropic theory, Politano et al., 1989, 
for the numerical confirmation of such a model in the two-dimensional isotropic case, and for the three- 
dimensional case, Politano et al., 1995). A small amount of correlation alters the energy spectrum only 
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slightly, making data analysis difficult. 

Among the other parameters of potential importance in the study of magnetohydrodynamic turbulence 
is the magnetic Prandtl number ly/r], the presence of sizable global (as opposed to local) correlation 
between the velocity and the magnetic field, the ratio of magnetic to kinetic energy, and the effect of 
an externally imposed uniform magnetic field. Compressibility, rotation, stratification are other external 
agents of imposed anisotropies, and important for plasma experiments, the presence of small-scale kinetic 
effects such as a Hall current or pressure anisotropies may also alter the dynamics of MHD flows. Using 
imposed symmetries to enhance the numerically attainable Reynolds number at a given cost may be a 
promising venue in many of these cases. 
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